PAINTstudies_AminoAcidStudies
Dr. Eva Caamaño-Gutiérrez - Computational Biology Facility
2023-12-21
We start by loading the packages and functions needed
library(limma)
library(samr)
library(psych)
library(gplots)
library(ggpubr)
library(conover.test)
library(gtools) #for stars.pval
source("PAINTfunctions.R")1. Input data and some exploratory analysis
Upon exploratory analysis using all the dataset we do not observe any particular structure.
aaPaint<-read.csv("data/PAINTStudiesBiochemDat.csv")
aaPaint$TimePoint<-factor(aaPaint$TimePoint,levels=c("Day3","Day10"),ordered = T)
str(aaPaint)## 'data.frame': 140 obs. of 31 variables:
## $ Study : chr "PAINT18" "PAINT18" "PAINT18" "PAINT18" ...
## $ Group : chr "Intervention" "Intervention" "Intervention" "Intervention" ...
## $ LevelIntervention: chr "High-18" "High-18" "High-18" "High-18" ...
## $ TimePoint : Ord.factor w/ 2 levels "Day3"<"Day10": 2 2 2 2 2 2 2 2 2 2 ...
## $ PatientID : int 1821 1822 1823 1824 1825 1826 1827 1828 1830 1851 ...
## $ Phe : int 87 86 67 96 45 51 65 53 76 52 ...
## $ Val : int 231 211 168 217 93 115 161 119 175 118 ...
## $ Leu : int 175 171 120 192 75 109 132 81 164 86 ...
## $ Iso : int 71 72 43 74 43 52 50 42 61 49 ...
## $ Lys : int 306 340 213 370 134 180 207 195 236 178 ...
## $ Met : int 34 50 27 35 17 10 19 27 26 34 ...
## $ Thr : int 190 296 434 621 66 309 207 194 290 56 ...
## $ His : int 90 120 114 127 62 63 71 75 97 70 ...
## $ Try : int 19 20 26 12 19 7 11 21 15 28 ...
## $ Total.EAA : int 1203 1366 1212 1744 554 896 923 807 1140 671 ...
## $ Tyr : int 23 49 85 45 83 13 23 110 24 92 ...
## $ Gln : int 499 724 566 743 370 255 319 365 452 342 ...
## $ Arg : int 146 253 154 182 30 117 100 54 162 61 ...
## $ Cys : int 17 30 11 14 24 7 11 29 17 20 ...
## $ Gly : int 321 474 354 399 229 147 272 284 328 280 ...
## $ Pro : int 467 471 285 607 175 262 184 200 252 197 ...
## $ Glu : int 132 125 61 167 84 96 62 62 151 150 ...
## $ Asn : int 21 28 28 38 32 15 23 52 21 47 ...
## $ Asp : int 45 47 17 55 30 28 23 21 40 38 ...
## $ Ala : int 396 595 418 885 241 214 210 221 305 307 ...
## $ Ser : int 223 285 282 250 146 138 149 151 201 164 ...
## $ Tau : int 88 385 116 74 155 40 133 242 153 369 ...
## $ Orn : int 416 666 297 522 52 289 320 97 396 82 ...
## $ Cit : int 11 11 6 6 14 9 12 20 20 11 ...
## $ Ammonia : int 27 74 54 9 67 21 79 67 78 79 ...
## $ ArgIVIntake : num 452.7 655.7 254.5 95.9 0 ...#PCA plots (not taking into account IVArg)
do_PCA_Plot_Nov16((aaPaint[,6:30]),groups = aaPaint$LevelIntervention,legendName = "Arg \n Suppl",scale=F)
do_PCA_Plot_Nov16((aaPaint[,6:30]),groups = aaPaint$TimePoint,legendName = "time",scale=F)
do_PCA_Plot_Nov16((aaPaint[,6:30]),groups = paste(aaPaint$LevelIntervention,aaPaint$TimePoint,sep=""),legendName = "",scale=T)
#PCA plot of day 10 coloured by intervention - no structure appreciated
do_PCA_Plot_Nov16((aaPaint[aaPaint$TimePoint=="Day10",6:30]),groups = aaPaint[aaPaint$TimePoint=="Day10","LevelIntervention"],legendName = "",scale=F)
2. Study correlation of plasma Arginine with IV Arginine intake
Using only the data available for patients on the intervention for PAINT18. We can see from the analysis shown below that while at day 3 there is no correlation between intake and levels of arginine, by day 10 there is a strong correlation (Spearman correlation 0.73 with p value <0.001). It is interesting to see that the patients that were on enteral feeds (zero levels at Day 10 for IV Arginine), show relatively low levels of Arginine. Note patient 1824 died the day after D10 sample.
ArgCorre<-aaPaint[!is.na(aaPaint$ArgIVIntake),]
ArgCorre$TimePoint<-factor(ArgCorre$TimePoint,levels=c("Day3","Day10"),ordered = T)
ggplot(data=ArgCorre,aes(x=Arg,y=ArgIVIntake))+
geom_point()+ggrepel::geom_label_repel(aes(label=PatientID),box.padding=0.35,point.padding = 0.5,segment.color="grey50")+facet_wrap(~TimePoint,scales="free_x")+
geom_smooth(aes(x=Arg,y=ArgIVIntake),method='lm')+
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 10))+
xlab("Arginine levels plasma")+ylab("Arginine IV Intake")+
stat_cor(aes(x=Arg,y=ArgIVIntake),method="spearman", col="blue",label.y = 775)+
theme_bw(base_size = 16)
3. Statistical analysis of total essential aminoacids, plasma arginine and amonia
We start by producing some visualisations
aaPaintM<-melt(aaPaint[,-ncol(aaPaint)],id.vars=c("Study","Group","LevelIntervention","TimePoint","PatientID"))
aaPaintM$variable<-as.character(aaPaintM$variable)
aaPaintM_EAA_Arg_Am<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia"),]
ggplot(data=aaPaintM_EAA_Arg_Am)+
geom_boxplot(aes(x=LevelIntervention, y=value,fill=TimePoint),col="black")+
geom_point(aes(x=LevelIntervention, y=value,fill=TimePoint,shape=Study),position=position_jitterdodge(dodge.width=0.9),col="black",alpha=0.4,pch=21)+
facet_wrap(~variable,scales="free")+
ylab("")+
scale_fill_manual(values = c("#DCE319FF","#33638DFF"))+
scale_colour_manual(values = c("#DCE319FF","#33638DFF"))+
scale_shape_manual(values=c(21,22,23))+
theme_bw(base_size = 16)+theme(legend.position = "bottom")
aaPaintM_sel<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Cit","Orn"),]
ggplot(data=aaPaintM_sel)+
geom_boxplot(aes(x=LevelIntervention, y=value,fill=TimePoint),col="black")+
geom_point(aes(x=LevelIntervention, y=value,fill=TimePoint,shape=Study),position=position_jitterdodge(dodge.width=0.9),col="black",alpha=0.4,pch=21)+
facet_wrap(~variable,scales="free")+
ylab("")+
scale_fill_manual(values = c("#DCE319FF","#33638DFF"))+
scale_colour_manual(values = c("#DCE319FF","#33638DFF"))+
scale_shape_manual(values=c(21,22,23))+
theme_bw(base_size = 16)+theme(legend.position = "bottom")
Comparisons for Day 10
Upon data exploration the statistical test chosen to look at the differences at day 10 for each metabolite is a Kruskal-Wallis, corrected for FDR with Benjamini-Hochberg followed by a Conover-Iman test as post-Hoc. Significant aminoacids shown in boxplots with standard significance values.
qqnorm(aaPaint$Arg)#no
shapiro.test(aaPaint$Arg)##
## Shapiro-Wilk normality test
##
## data: aaPaint$Arg
## W = 0.7976, p-value = 1.255e-12qqnorm(aaPaint$Ammonia)#yes
shapiro.test(aaPaint$Ammonia)##
## Shapiro-Wilk normality test
##
## data: aaPaint$Ammonia
## W = 0.98934, p-value = 0.362qqnorm(aaPaint$Total.EAA)#No
shapiro.test(aaPaint$Total.EAA)##
## Shapiro-Wilk normality test
##
## data: aaPaint$Total.EAA
## W = 0.9713, p-value = 0.004793#with the above prelim tests it is better to undertake a non-parametric test#we'll start by doing a Kruskall Wallis at all the day 10 measurements followed by post-hoc with conover test
aaPaint_D10<-aaPaint[aaPaint$TimePoint=="Day10",]
KW_all_pvals<-apply(aaPaint_D10[,6:30],2,function(x){kruskal.test(x,aaPaint_D10$LevelIntervention)$p.value})
#we correct for FDR and we select the significant for the posthoc
KW_all_pvals_adj<-p.adjust(KW_all_pvals,method = "BH")
namesSig<-names(KW_all_pvals_adj[KW_all_pvals_adj<0.05])#we get 8 significant
forConover<-aaPaint_D10[,colnames(aaPaint_D10)%in%namesSig]
KW_conoverIman_postHoc<-do.call(rbind,lapply(lapply(forConover,conover.test, g=aaPaint_D10$LevelIntervention,list=T,method="BH",label=T,table=F),as.data.frame))## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 18.9695, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -2.194453 (0.0159)*
## High-18 - Standard_6 : -4.970219 (0.0000)*
## Moderate - Standard_6 : -3.134436 (0.0019)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 10.32, df = 2, p-value = 0.01
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -0.649611 (0.2591)
## High-18 - Standard_6 : -3.092535 (0.0044)*
## Moderate - Standard_6 : -2.747346 (0.0058)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 7.3808, df = 2, p-value = 0.02
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : 0.064360 (0.4744)
## High-18 - Standard_6 : -2.230323 (0.0219)*
## Moderate - Standard_6 : -2.574716 (0.0185)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 11.8363, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -3.575898 (0.0010)*
## High-18 - Standard_6 : -1.283648 (0.1019)
## Moderate - Standard_6 : 2.541188 (0.0101)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 8.5445, df = 2, p-value = 0.01
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -3.059427 (0.0048)*
## High-18 - Standard_6 : -2.135894 (0.0273)
## Moderate - Standard_6 : 1.009632 (0.1582)
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 20.9383, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : 1.772715 (0.0405)
## High-18 - Standard_6 : 5.174610 (0.0000)*
## Moderate - Standard_6 : 3.833430 (0.0002)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 19.2275, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -3.725538 (0.0003)*
## High-18 - Standard_6 : 0.568754 (0.2857)
## Moderate - Standard_6 : 4.786731 (0.0000)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 12.9503, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -3.933887 (0.0003)*
## High-18 - Standard_6 : -2.108928 (0.0291)*
## Moderate - Standard_6 : 2.013617 (0.0241)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 7.5947, df = 2, p-value = 0.02
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -1.124172 (0.1325)
## High-18 - Standard_6 : -2.798143 (0.0101)*
## Moderate - Standard_6 : -1.888525 (0.0476)
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 19.5445, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -2.005475 (0.0246)*
## High-18 - Standard_6 : 2.696496 (0.0067)*
## Moderate - Standard_6 : 5.198173 (0.0000)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 15.5167, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -0.364130 (0.3585)
## High-18 - Standard_6 : -3.747581 (0.0003)*
## Moderate - Standard_6 : -3.800377 (0.0005)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 9.2122, df = 2, p-value = 0.01
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -3.153203 (0.0037)*
## High-18 - Standard_6 : -1.373626 (0.0871)
## Moderate - Standard_6 : 1.969532 (0.0399)
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 20.0594, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -0.935407 (0.1765)
## High-18 - Standard_6 : 3.623966 (0.0004)*
## Moderate - Standard_6 : 5.067460 (0.0000)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 16.4153, df = 2, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -0.034309 (0.4864)
## High-18 - Standard_6 : 3.653762 (0.0004)*
## Moderate - Standard_6 : 4.138753 (0.0002)*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 7.7156, df = 2, p-value = 0.02
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -------------------------------------------
## High-18 - Moderate : -2.759684 (0.0113)*
## High-18 - Standard_6 : -2.381707 (0.0152)*
## Moderate - Standard_6 : 0.441281 (0.3303)
##
## alpha = 0.05
## Reject Ho if p <= alpha/2KW_conoverIman_postHoc$AminoAcid<-strsplit2(rownames(KW_conoverIman_postHoc),"[.]")[,1]This is the table with the post-hoc analysis for the significant aminoacids.
KW_conoverIman_postHoc[KW_conoverIman_postHoc$AminoAcid%in%c("Total"),"AminoAcid"]<-"Total.EAA"
KW_conoverIman_postHoc[,c(6,5,2,4)]## AminoAcid comparisons T P.adjusted
## Phe.1 Phe High-18 - Moderate -2.19445336 1.589113e-02
## Phe.2 Phe High-18 - Standard_6 -4.97021968 7.721494e-06
## Phe.3 Phe Moderate - Standard_6 -3.13443624 1.937981e-03
## Val.1 Val High-18 - Moderate -0.64961127 2.591162e-01
## Val.2 Val High-18 - Standard_6 -3.09253593 4.384765e-03
## Val.3 Val Moderate - Standard_6 -2.74734601 5.824426e-03
## Leu.1 Leu High-18 - Moderate 0.06436025 4.744405e-01
## Leu.2 Leu High-18 - Standard_6 -2.23032310 2.188851e-02
## Leu.3 Leu Moderate - Standard_6 -2.57471604 1.847714e-02
## Thr.1 Thr High-18 - Moderate -3.57589808 9.972381e-04
## Thr.2 Thr High-18 - Standard_6 -1.28364837 1.019112e-01
## Thr.3 Thr Moderate - Standard_6 2.54118872 1.008235e-02
## Total.EAA.1 Total.EAA High-18 - Moderate -3.05942756 4.830209e-03
## Total.EAA.2 Total.EAA High-18 - Standard_6 -2.13589452 2.734474e-02
## Total.EAA.3 Total.EAA Moderate - Standard_6 1.00963286 1.582075e-01
## Arg.1 Arg High-18 - Moderate 1.77271522 4.048072e-02
## Arg.2 Arg High-18 - Standard_6 5.17461063 3.578830e-06
## Arg.3 Arg Moderate - Standard_6 3.83343065 2.157958e-04
## Cys.1 Cys High-18 - Moderate -3.72553847 3.076751e-04
## Cys.2 Cys High-18 - Standard_6 0.56875407 2.857415e-01
## Cys.3 Cys Moderate - Standard_6 4.78673184 1.524334e-05
## Gly.1 Gly High-18 - Moderate -3.93388718 3.087180e-04
## Gly.2 Gly High-18 - Standard_6 -2.10892838 2.910623e-02
## Gly.3 Gly Moderate - Standard_6 2.01361725 2.409668e-02
## Pro.1 Pro High-18 - Moderate -1.12417239 1.325372e-01
## Pro.2 Pro High-18 - Standard_6 -2.79814379 1.013419e-02
## Pro.3 Pro Moderate - Standard_6 -1.88852595 4.756429e-02
## Asn.1 Asn High-18 - Moderate -2.00547530 2.460601e-02
## Asn.2 Asn High-18 - Standard_6 2.69649638 6.733013e-03
## Asn.3 Asn Moderate - Standard_6 5.19817339 3.474889e-06
## Asp.1 Asp High-18 - Moderate -0.36413094 3.584710e-01
## Asp.2 Asp High-18 - Standard_6 -3.74758181 2.862929e-04
## Asp.3 Asp Moderate - Standard_6 -3.80037745 4.814110e-04
## Ser.1 Ser High-18 - Moderate -3.15320345 3.666462e-03
## Ser.2 Ser High-18 - Standard_6 -1.37362619 8.713847e-02
## Ser.3 Ser Moderate - Standard_6 1.96953297 3.986740e-02
## Tau.1 Tau High-18 - Moderate -0.93540725 1.765473e-01
## Tau.2 Tau High-18 - Standard_6 3.62396628 4.315131e-04
## Tau.3 Tau Moderate - Standard_6 5.06746037 5.512059e-06
## Orn.1 Orn High-18 - Moderate -0.03430907 4.863679e-01
## Orn.2 Orn High-18 - Standard_6 3.65376272 3.883769e-04
## Orn.3 Orn Moderate - Standard_6 4.13875353 1.537773e-04
## Cit.1 Cit High-18 - Moderate -2.75968448 1.134998e-02
## Cit.2 Cit High-18 - Standard_6 -2.38170764 1.519900e-02
## Cit.3 Cit Moderate - Standard_6 0.44128108 3.302602e-01write.csv(KW_conoverIman_postHoc,"D10_KW_ConoverImanPH.csv")Now we generate a plot with the results of the significant aminoacids.
aaPaintM_D10_sig<-aaPaintM[aaPaintM$TimePoint=="Day10" & aaPaintM$variable%in%namesSig,]
colnames(aaPaintM_D10_sig)[6]<-"aminoacid"
KW_conoverIman_postHoc$group1<-strsplit2(KW_conoverIman_postHoc$comparisons," - ")[,1]
KW_conoverIman_postHoc$group2<-strsplit2(KW_conoverIman_postHoc$comparisons," - ")[,2]
ypos<-apply(aaPaint[,namesSig],2,max,na.rm=T)*1.12
KW_conoverIman_postHoc$y.position<-ypos[match(KW_conoverIman_postHoc$AminoAcid,names(ypos))]
KW_conoverIman_postHoc$StatSig<-stars.pval(KW_conoverIman_postHoc$P.adjusted)
KW_conoverIman_postHoc$StatSig[KW_conoverIman_postHoc$StatSig==" "]<-"NS"
KW_conoverIman_postHoc[,c(6,5,2,4)]## AminoAcid comparisons T P.adjusted
## Phe.1 Phe High-18 - Moderate -2.19445336 1.589113e-02
## Phe.2 Phe High-18 - Standard_6 -4.97021968 7.721494e-06
## Phe.3 Phe Moderate - Standard_6 -3.13443624 1.937981e-03
## Val.1 Val High-18 - Moderate -0.64961127 2.591162e-01
## Val.2 Val High-18 - Standard_6 -3.09253593 4.384765e-03
## Val.3 Val Moderate - Standard_6 -2.74734601 5.824426e-03
## Leu.1 Leu High-18 - Moderate 0.06436025 4.744405e-01
## Leu.2 Leu High-18 - Standard_6 -2.23032310 2.188851e-02
## Leu.3 Leu Moderate - Standard_6 -2.57471604 1.847714e-02
## Thr.1 Thr High-18 - Moderate -3.57589808 9.972381e-04
## Thr.2 Thr High-18 - Standard_6 -1.28364837 1.019112e-01
## Thr.3 Thr Moderate - Standard_6 2.54118872 1.008235e-02
## Total.EAA.1 Total.EAA High-18 - Moderate -3.05942756 4.830209e-03
## Total.EAA.2 Total.EAA High-18 - Standard_6 -2.13589452 2.734474e-02
## Total.EAA.3 Total.EAA Moderate - Standard_6 1.00963286 1.582075e-01
## Arg.1 Arg High-18 - Moderate 1.77271522 4.048072e-02
## Arg.2 Arg High-18 - Standard_6 5.17461063 3.578830e-06
## Arg.3 Arg Moderate - Standard_6 3.83343065 2.157958e-04
## Cys.1 Cys High-18 - Moderate -3.72553847 3.076751e-04
## Cys.2 Cys High-18 - Standard_6 0.56875407 2.857415e-01
## Cys.3 Cys Moderate - Standard_6 4.78673184 1.524334e-05
## Gly.1 Gly High-18 - Moderate -3.93388718 3.087180e-04
## Gly.2 Gly High-18 - Standard_6 -2.10892838 2.910623e-02
## Gly.3 Gly Moderate - Standard_6 2.01361725 2.409668e-02
## Pro.1 Pro High-18 - Moderate -1.12417239 1.325372e-01
## Pro.2 Pro High-18 - Standard_6 -2.79814379 1.013419e-02
## Pro.3 Pro Moderate - Standard_6 -1.88852595 4.756429e-02
## Asn.1 Asn High-18 - Moderate -2.00547530 2.460601e-02
## Asn.2 Asn High-18 - Standard_6 2.69649638 6.733013e-03
## Asn.3 Asn Moderate - Standard_6 5.19817339 3.474889e-06
## Asp.1 Asp High-18 - Moderate -0.36413094 3.584710e-01
## Asp.2 Asp High-18 - Standard_6 -3.74758181 2.862929e-04
## Asp.3 Asp Moderate - Standard_6 -3.80037745 4.814110e-04
## Ser.1 Ser High-18 - Moderate -3.15320345 3.666462e-03
## Ser.2 Ser High-18 - Standard_6 -1.37362619 8.713847e-02
## Ser.3 Ser Moderate - Standard_6 1.96953297 3.986740e-02
## Tau.1 Tau High-18 - Moderate -0.93540725 1.765473e-01
## Tau.2 Tau High-18 - Standard_6 3.62396628 4.315131e-04
## Tau.3 Tau Moderate - Standard_6 5.06746037 5.512059e-06
## Orn.1 Orn High-18 - Moderate -0.03430907 4.863679e-01
## Orn.2 Orn High-18 - Standard_6 3.65376272 3.883769e-04
## Orn.3 Orn Moderate - Standard_6 4.13875353 1.537773e-04
## Cit.1 Cit High-18 - Moderate -2.75968448 1.134998e-02
## Cit.2 Cit High-18 - Standard_6 -2.38170764 1.519900e-02
## Cit.3 Cit Moderate - Standard_6 0.44128108 3.302602e-01write.csv(KW_conoverIman_postHoc[,c(6,5,2,4)],"ConoverImanResults.csv")
for(i in namesSig){
paa<-ggplot(data=aaPaintM_D10_sig[aaPaintM_D10_sig$aminoacid==i,],aes(x=LevelIntervention))+
geom_boxplot(aes(x=LevelIntervention, y=value),col="black")+
geom_jitter(aes(x=LevelIntervention, y=value),col="black",alpha=0.6,pch=16)+
facet_wrap(~aminoacid,scales="free")+
ylab("")+xlab("")+
theme_classic(base_size = 16)+theme(legend.position = "bottom")
paa2<-paa+stat_pvalue_manual(data=KW_conoverIman_postHoc[KW_conoverIman_postHoc$AminoAcid==i,],label="StatSig",step.increase = c(0,0.1,0.1))
ggsave(plot = paa2,filename = paste("FigBiochem/",i,"_Day10_boxplot.png",sep = ""))
print(paa2)
}
Comparisons for day 3
While the first strategy included looking at comparisons between the three groups, as clarified by the clinical team not all the studies are equivalent as some started the intervention after the sample was taken (PAINT 18). Hence Day 3 comparisons will be done in a 3-side non parametric test
Day 3 PAINT 18 data
#for PAINT18
aaPaint_D3_PAINT18<-aaPaint[aaPaint$TimePoint=="Day3" & aaPaint$Study=="PAINT18",]
D3PAINT18_stats_pvals<-sapply(6:30,function(x){
wilcox.test(x = aaPaint_D3_PAINT18[aaPaint_D3_PAINT18$LevelIntervention=="Standard_6",x],
y = aaPaint_D3_PAINT18[aaPaint_D3_PAINT18$LevelIntervention=="High-18",x])$p.value})
D3PAINT18_stats<-data.frame("aminoacid"=colnames(aaPaint_D3_PAINT18)[6:30],"pvalue"=D3PAINT18_stats_pvals,"Padj"=p.adjust(D3PAINT18_stats_pvals,method = "BH"))
write.csv(D3PAINT18_stats,"D3PAINT18_CvsIComparison_SuppTable3.csv")Day 3 PAINT data
aaPaint_D3_PAINT<-aaPaint[aaPaint$TimePoint=="Day3" & aaPaint$Study%in%c("PAINT","PAINT-NH3"),]
D3PAINT_stats_pvals<-sapply(6:30,function(x){
wilcox.test(x = aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",x],
y = aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",x])$p.value})
D3PAINT_stats<-data.frame("aminoacid"=colnames(aaPaint_D3_PAINT)[6:30],"pvalue"=D3PAINT_stats_pvals,"Padj"=p.adjust(D3PAINT_stats_pvals,method = "BH"))
D3PAINT_stats$medianControl<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",6:30],2,median,na.rm=T)
D3PAINT_stats$medianInterv<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",6:30],2,median,na.rm=T)
D3PAINT_stats$IQR_Control<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",6:30],2,IQR,na.rm=T)
D3PAINT_stats$IQR_Interv<-apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",6:30],2,IQR,na.rm=T)
write.csv(D3PAINT_stats,"D3PAINT_CvsIComparison_SuppTable2.csv")
apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Standard_6",6:30],2,quantile,na.rm=T)## Phe Val Leu Iso Lys Met Thr His Try Total.EAA Tyr
## 0% 59.00 113.00 87.00 13.00 135.00 13.00 135.00 54.0 9.0 788.00 11.00
## 25% 66.25 160.25 131.25 46.00 172.75 22.75 260.25 76.5 14.0 1030.50 32.25
## 50% 88.50 180.00 145.00 50.50 212.50 27.50 325.00 85.5 24.5 1220.00 79.00
## 75% 98.00 204.25 157.00 64.25 275.25 40.25 511.00 106.0 32.5 1427.75 143.25
## 100% 124.00 272.00 204.00 80.00 424.00 50.00 626.00 141.0 41.0 1765.00 657.00
## Gln Arg Cys Gly Pro Glu Asn Asp Ala Ser Tau
## 0% 157.00 5.00 5.00 237.0 215.00 60.00 15.00 135.00 197.0 90.00 17.0
## 25% 358.75 11.75 10.75 275.5 321.50 101.25 19.75 172.75 275.5 166.50 42.5
## 50% 442.50 19.50 16.50 351.0 339.00 121.00 24.00 212.50 389.5 225.00 57.0
## 75% 527.00 36.25 22.75 447.5 440.25 177.75 33.25 275.25 428.5 280.75 86.0
## 100% 664.00 61.00 42.00 748.0 492.00 267.00 47.00 424.00 556.0 338.00 187.0
## Orn Cit Ammonia
## 0% 22.00 11.00 28.00
## 25% 68.25 12.75 43.50
## 50% 91.00 14.50 52.50
## 75% 120.00 19.50 63.25
## 100% 184.00 27.00 78.00apply(aaPaint_D3_PAINT[aaPaint_D3_PAINT$LevelIntervention=="Moderate",6:30],2,quantile,na.rm=T)## Phe Val Leu Iso Lys Met Thr His Try Total.EAA Tyr
## 0% 63 109.00 41.00 30.00 136.00 13.00 113.00 55.0 11.0 696.0 17.0
## 25% 74 137.75 115.75 43.75 173.50 22.75 208.75 68.0 16.5 929.0 40.5
## 50% 85 179.50 145.50 51.50 206.00 28.00 267.50 90.0 24.0 1072.0 74.0
## 75% 104 201.25 155.50 58.00 250.75 35.00 293.75 108.5 32.5 1290.5 105.5
## 100% 129 331.00 289.00 103.00 365.00 57.00 463.00 137.0 56.0 1823.0 332.0
## Gln Arg Cys Gly Pro Glu Asn Asp Ala Ser Tau Orn
## 0% 214.00 13.00 3.0 209.0 181.0 74.0 14.0 19.00 183.0 116.00 13.00 60.0
## 25% 348.75 28.50 7.0 274.5 243.5 107.5 21.5 25.75 276.5 153.50 40.25 112.5
## 50% 467.50 40.00 14.5 327.5 276.0 137.5 27.0 30.50 320.0 213.50 77.50 165.5
## 75% 540.25 57.25 23.5 411.0 344.5 171.5 33.0 39.25 459.0 269.75 120.50 227.0
## 100% 676.00 164.00 49.0 510.0 613.0 273.0 78.0 72.00 656.0 390.00 264.00 468.0
## Cit Ammonia
## 0% 8.00 9.00
## 25% 13.00 38.25
## 50% 14.00 41.00
## 75% 17.25 46.00
## 100% 26.00 64.00Day 30 PAINT 18 data
#for PAINT18
aaPaintWithD30<-read.csv("data/PAINTStudiesBiochemDat_withD30.csv")
aaPaint_D30_PAINT18<-aaPaintWithD30[aaPaintWithD30$TimePoint=="Day30" & aaPaintWithD30$Study=="PAINT18",]
D30PAINT18_stats_pvals<-sapply(6:30,function(x){
wilcox.test(x = aaPaint_D30_PAINT18[aaPaint_D30_PAINT18$LevelIntervention=="Standard_6",x],
y = aaPaint_D30_PAINT18[aaPaint_D30_PAINT18$LevelIntervention=="High-18",x],na.rm=T)$p.value})
D30PAINT18_stats<-data.frame("aminoacid"=colnames(aaPaint_D30_PAINT18)[6:30],"pvalue"=D30PAINT18_stats_pvals,"Padj"=p.adjust(D30PAINT18_stats_pvals,method = "BH"))
write.csv(D30PAINT18_stats,"D30PAINT18_CvsIComparison_SuppTable4.csv")Below find preparation of some of the figures asked by the medical team for the manuscript:
#boxplots comparing key variables at day 3
aaPaintM<-melt(aaPaint[,-ncol(aaPaint)],id.vars=c("Study","Group","LevelIntervention","TimePoint","PatientID"))
aaPaintM$variable<-as.character(aaPaintM$variable)
aaPaintM_EAA_Arg_Am_Day3<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Gln","Cit","Orn") & aaPaintM$TimePoint=="Day3",]
boxplotD3<-ggplot(data=aaPaintM_EAA_Arg_Am_Day3)+
geom_boxplot(aes(x=LevelIntervention, y=value),col="black",outlier.color = "white")+
geom_point(aes(x=LevelIntervention, y=value,shape=Study),alpha=0.6)+
facet_wrap(~variable,scales="free")+
ylab("Abundance")+ xlab("Intervention")+
theme_classic(base_size = 12)+theme(legend.position = "bottom")
boxplotD3
#Day 10
aaPaintM_EAA_Arg_Am_Day10<-aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Gln","Cit","Orn") & aaPaintM$TimePoint=="Day10",]
aaPaintM_EAA_Arg_Am_Day10$variable<-factor(aaPaintM_EAA_Arg_Am_Day10$variable,ordered=T,levels=c("Arg","Total.EAA","Ammonia","Orn","Cit","Gln"))
aaPaintM_EAA_Arg_Am_Day10$Study[aaPaintM_EAA_Arg_Am_Day10$Study%in%c("PAINT","PAINT-NH3")]<-"PAINT"
boxplotD10<-ggplot(data=aaPaintM_EAA_Arg_Am_Day10)+
geom_boxplot(aes(x=LevelIntervention, y=value),col="black",outlier.color = "white")+
geom_point(aes(x=LevelIntervention, y=value,shape=Study),alpha=0.6)+
facet_wrap(~variable,scales="free")+
ylab("Abundance")+ xlab("Intervention")+
theme_classic(base_size = 12)+theme(legend.position = "bottom")
boxplotD10
ggsave(plot=boxplotD10,filename = "D10_SelectedAminoAcids.png",device = "png",dpi = 300,width=8,units = "in")
ggplot(data=aaPaintM[aaPaintM$variable%in%c("Total.EAA","Arg","Ammonia","Gln","Cit","Orn"),])+
geom_boxplot(aes(x=LevelIntervention, y=value, fill=TimePoint),col="black",outlier.color = "white")+
geom_point(aes(x=LevelIntervention, y=value,shape=Study,col=TimePoint,fill=TimePoint),alpha=0.6,position=position_jitterdodge(dodge.width=0.9))+
scale_fill_manual(values = c("darkgrey","skyblue"))+
scale_colour_manual(values = c("darkgrey","skyblue"))+
facet_wrap(~variable,scales="free")+
ylab("Abundance")+ xlab("Intervention")+
theme_classic(base_size = 12)+theme(legend.position = "bottom")